home *** CD-ROM | disk | FTP | other *** search
/ io Programmo 60 / IOPROG_60.ISO / soft / c++ / gsl-1.1.1-setup.exe / {app} / src / specfunc / clausen.c < prev    next >
Encoding:
C/C++ Source or Header  |  2002-04-18  |  2.8 KB  |  112 lines
  1. /* specfunc/clausen.c
  2.  * 
  3.  * Copyright (C) 1996, 1997, 1998, 1999, 2000 Gerard Jungman
  4.  * 
  5.  * This program is free software; you can redistribute it and/or modify
  6.  * it under the terms of the GNU General Public License as published by
  7.  * the Free Software Foundation; either version 2 of the License, or (at
  8.  * your option) any later version.
  9.  * 
  10.  * This program is distributed in the hope that it will be useful, but
  11.  * WITHOUT ANY WARRANTY; without even the implied warranty of
  12.  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13.  * General Public License for more details.
  14.  * 
  15.  * You should have received a copy of the GNU General Public License
  16.  * along with this program; if not, write to the Free Software
  17.  * Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
  18.  */
  19.  
  20. /* Author:  G. Jungman */
  21.  
  22. #include <config.h>
  23. #include <gsl/gsl_math.h>
  24. #include <gsl/gsl_errno.h>
  25. #include <gsl/gsl_sf_trig.h>
  26. #include <gsl/gsl_sf_clausen.h>
  27.  
  28. #include "chebyshev.h"
  29. #include "cheb_eval.c"
  30.  
  31. static double aclaus_data[15] = {
  32.   2.142694363766688447e+00,
  33.   0.723324281221257925e-01,
  34.   0.101642475021151164e-02,
  35.   0.3245250328531645e-04,
  36.   0.133315187571472e-05,
  37.   0.6213240591653e-07,
  38.   0.313004135337e-08,
  39.   0.16635723056e-09,
  40.   0.919659293e-11,
  41.   0.52400462e-12,
  42.   0.3058040e-13,
  43.   0.18197e-14,
  44.   0.1100e-15,
  45.   0.68e-17,
  46.   0.4e-18
  47. };
  48. static cheb_series aclaus_cs = {
  49.   aclaus_data,
  50.   14,
  51.   -1, 1,
  52.   8  /* FIXME:  this is a guess, correct value needed here BJG */
  53. };
  54.  
  55.  
  56. /*-*-*-*-*-*-*-*-*-*-*-* Functions with Error Codes *-*-*-*-*-*-*-*-*-*-*-*/
  57.  
  58. int gsl_sf_clausen_e(double x, gsl_sf_result *result)
  59. {
  60.   const double x_cut = M_PI * GSL_SQRT_DBL_EPSILON;
  61.  
  62.   double sgn = 1.0;
  63.   int status_red;
  64.  
  65.   if(x < 0.0) {
  66.     x   = -x;
  67.     sgn = -1.0;
  68.   }
  69.  
  70.   /* Argument reduction to [0, 2pi) */
  71.   status_red = gsl_sf_angle_restrict_pos_e(&x);
  72.  
  73.   /* Further reduction to [0,pi) */
  74.   if(x > M_PI) {
  75.     /* simulated extra precision: 2PI = p0 + p1 */
  76.     const double p0 = 6.28125;
  77.     const double p1 = 0.19353071795864769253e-02;
  78.     x = (p0 - x) + p1;
  79.     sgn = -sgn;
  80.   }
  81.  
  82.   if(x == 0.0) {
  83.     result->val = 0.0;
  84.     result->err = 0.0;
  85.   }
  86.   else if(x < x_cut) {
  87.     result->val = x * (1.0 - log(x));
  88.     result->err = x * GSL_DBL_EPSILON;
  89.   }
  90.   else {
  91.     const double t = 2.0*(x*x / (M_PI*M_PI) - 0.5);
  92.     gsl_sf_result result_c;
  93.     cheb_eval_e(&aclaus_cs, t, &result_c);
  94.     result->val = x * (result_c.val - log(x));
  95.     result->err = x * (result_c.err + GSL_DBL_EPSILON);
  96.   }
  97.  
  98.   result->val *= sgn;
  99.  
  100.   return status_red;
  101. }
  102.  
  103.  
  104. /*-*-*-*-*-*-*-*-*-* Functions w/ Natural Prototypes *-*-*-*-*-*-*-*-*-*-*/
  105.  
  106. #include "eval.h"
  107.  
  108. double gsl_sf_clausen(const double x)
  109. {
  110.   EVAL_RESULT(gsl_sf_clausen_e(x, &result));
  111. }
  112.